Names: Ricardo Leão Cardoso, Annabelle Ducret, Alexandrine Diff

PCA for Source Separation of Abdominal ECG Signals¶

Introduction¶

In this exercise we use PCA for the separation of maternal and fetal electrocardiography (ECG) signals in abdominal ECG (aECG) data recorded on the belly of a pregnant woman. Due to the low signal strengh of fetal ECG (fECG) signals it is an "algorithmic challenge" to properly separate fECG from much stronger maternal ECG (mECG) signals [1].

The present example uses a simplified version of the method proposed by Varanini et al. [2].

References¶

[1] R. Kahankova et al., “A Review of Signal Processing Techniques for Non-Invasive Fetal Electrocardiography,” IEEE Reviews in Biomedical Engineering, vol. 13, pp. 51–73, 2020, doi: 10.1109/RBME.2019.2938061.

[2] M. Varanini, G. Tartarisco, L. Billeci, A. Macerata, G. Pioggia, and R. Balocchi, “An efficient unsupervised fetal QRS complex detection from abdominal maternal ECG,” Physiol. Meas., vol. 35, no. 8, pp. 1607–1619, Aug. 2014, doi: 10.1088/0967-3334/35/8/1607.

[3] Source of data: https://physionet.org/content/challenge-2013/1.0.0/

In [1]:
%matplotlib widget

import numpy as np
import pandas as pd
from scipy.signal import filtfilt, butter
from sklearn.decomposition import PCA

import plotly.graph_objects as go
from plotly.offline import init_notebook_mode, iplot
from plotly.subplots import make_subplots
init_notebook_mode(connected=True)  # initiate notebook for offline plot

from mqrs_utils import cancel_mqrs
from ecgdetectors import Detectors
In [3]:
# load abdominal ECG (aECG) data
# transformed from initial source: https://physionet.org/content/challenge-2013/1.0.0/set-a/a13.dat
filename = 'aecg_a13.hdf5'
aecg = pd.read_hdf(filename, key='signals').values
fs = 1000
t = np.arange(aecg.shape[0]) / fs
In [5]:
# bandpass filter data
b, a = butter(4, np.asarray([3, 45]), fs=fs, btype='bandpass')
aecg = filtfilt(b, a, aecg, axis=0)

# plot
fig = go.Figure()
for i in range(aecg.shape[1]):
    fig.add_trace(go.Scatter(x=t, y=aecg[:, i], name='AECG{:d}'.format(i)))
fig.update_xaxes(title='Time (s)')
fig.update_yaxes(title='AECG Amplitude (A.U.)')
fig.update_layout(title='Bandpass-Filtered AECG Signals')
fig.show()
In [7]:
# apply first PCA for enhancing maternal ECG component
pca1 = PCA()
pc1 = pca1.fit_transform(aecg)

# maternal ECG as the first principal component, note that this 
# remains a guess and would need to be automated in the final solution
maternal_ecg = pc1[:, 0]

# detect maternal QRS peaks
mqrs_peaks = Detectors(fs).engzee_detector(maternal_ecg)

# plot
fig = go.Figure()
for i in range(pc1.shape[1]):
    fig.add_trace(go.Scatter(x=t, y=pc1[:, i], name='PC1[:,{:d}]'.format(i)))
    if i == 0:
        fig.add_trace(go.Scatter(x=t[mqrs_peaks], y=pc1[mqrs_peaks, i], name='mQRS-Peaks', 
                                 mode='markers', marker_color='red', marker_symbol='circle-open'))
fig.update_xaxes(title='Time (s)')
fig.update_yaxes(title='PC1 (A.U.)')
fig.update_layout(title='Principal Components of First PCA Used to Enhance mECG Signal')
fig.show()
In [9]:
# remove maternal QRS complexes from signal to obtain a best possible fetal ECG signal
x_residual, mecg_estimations = cancel_mqrs(fs, pc1, np.asarray(mqrs_peaks))

# plot
fig = make_subplots(rows=3, cols=1, shared_xaxes=True)
fig.add_trace(go.Scatter(x=t, y=pc1[:,0], name='Maternal ECG'), row=1, col=1)
fig.add_trace(go.Scatter(x=t[mqrs_peaks], y=pc1[mqrs_peaks, 0], name='mQRS-Peaks', marker_color='red',
                         legendgroup='mQRS', mode='markers', marker_symbol='circle-open'), row=1, col=1)    
fig.add_trace(go.Scatter(x=t, y=mecg_estimations[:, 0], name='Interpolated mQRS Signal'), row=2, col=1)
fig.add_trace(go.Scatter(x=t[mqrs_peaks], y=mecg_estimations[mqrs_peaks, 0], name='mQRS-Peaks', marker_color='red', 
                         legendgroup='mQRS', showlegend=False, mode='markers', marker_symbol='circle-open'), row=2, col=1)    
fig.add_trace(go.Scatter(x=t, y=x_residual[:, 0], name='mQRS-free Signal'), row=3, col=1)
fig.add_trace(go.Scatter(x=t[mqrs_peaks], y=x_residual[mqrs_peaks, 0], name='mQRS-Peaks', marker_color='red',
                         legendgroup='mQRS', showlegend=False, mode='markers', marker_symbol='circle-open'), row=3, col=1)    
fig.update_xaxes(title='Time (s)', row=3, col=1)
fig.update_layout(title='Maternal QRS Cancellation')
fig.show()
In [11]:
# apply second PCA for enhancing fetal ECG component in residual signal
pca2 = PCA()
pc2 = pca2.fit_transform(x_residual)

# fetal ECG as the first principal component, note that this 
# remains a guess and would need to be automated in the final solution
fetal_ecg = pc2[:, 0]
# detect fetal QRS peaks
fqrs_peaks = Detectors(fs).engzee_detector(fetal_ecg)

# plot
fig = go.Figure()
for i in range(pc2.shape[1]):
    fig.add_trace(go.Scatter(x=t, y=pc2[:, i], name='PC2[:,{:d}]'.format(i)))
    if i == 0:
        fig.add_trace(go.Scatter(x=t[fqrs_peaks], y=pc2[fqrs_peaks, i], name='fQRS-Peaks', 
                                 mode='markers', marker_color='black', marker_symbol='triangle-down-open'))
fig.update_xaxes(title='Time (s)')
fig.update_yaxes(title='PC2 (A.U.)')
fig.update_layout(title='Principal Components of Second PCA Used to Enhance fECG Signal')
fig.show()
In [13]:
# plot for summarizing all
fig = make_subplots(rows=2, cols=1, shared_xaxes=True)
# maternal ECG with mQRS
fig.add_trace(go.Scatter(x=t, y=maternal_ecg, name='Maternal ECG'), row=1, col=1)
fig.add_trace(go.Scatter(x=t[mqrs_peaks], y=maternal_ecg[mqrs_peaks], name='mQRS-Peaks', 
                         marker_color='red', mode='markers', marker_symbol='circle-open'), row=1, col=1)   
# fetal ECG with fQRS
fig.add_trace(go.Scatter(x=t, y=fetal_ecg, name='Fetal ECG'), row=2, col=1)
fig.add_trace(go.Scatter(x=t[fqrs_peaks], y=fetal_ecg[fqrs_peaks], name='fQRS-Peaks', 
                         marker_color='black', mode='markers', marker_symbol='triangle-down-open'), row=2, col=1)   
fig.update_xaxes(title='Time (s)', row=2, col=1)
fig.update_layout(title='Maternal vs. Fetal ECG')
fig.show()

Exercise Questions¶

Please provide your answers directly below each question.

Question 1¶

Determine the maternal heart rate, both expressed in Hz and beats/min.

In [27]:
f=1/np.mean(np.diff(t[mqrs_peaks]))
fm= 1/np.median(np.diff(t[mqrs_peaks]))

print('mean of the maternal heart rate: ' f'Hz:{f}, bpm:{f*60}')
print('median of the maternal heart rate: ' f'Hz:{fm}, bpm:{fm*60}')

mean of the maternal heart rate: Hz:1.3593882752761255, bpm:81.56329651656753
median of the maternal heart rate: Hz:1.3210039630118866, bpm:79.2602377807132

Question 2¶

Determine the fetal heart rate, both expressed in Hz and beats/min.

In [29]:
f=1/np.mean(np.diff(t[fqrs_peaks]))
fm= 1/np.median(np.diff(t[fqrs_peaks]))

print('mean of the fetal heart rate: 'f'Hz:{f}, bpm:{f*60}')
print('median of the fetal heart rate: 'f'Hz:{fm}, bpm:{fm*60}')

mean of the fetal heart rate: Hz:1.7426613653667202, bpm:104.55968192200321
median of the fetal heart rate: Hz:2.0964360587001925, bpm:125.78616352201155

Question 3¶

Determine the follwing three values:

  • i) the average amplitude of the maternal QRS peaks (mQRS);
  • ii) the average amplitude of the fetal QRS peaks (fQRS);
  • iii) the ratio between the average amplitudes of i) mQRS and ii) fQRS peaks.
In [39]:
AmQRS=np.mean(maternal_ecg[mqrs_peaks])
AfQRS=np.mean(fetal_ecg[fqrs_peaks])
print(f'average amplitude of mQRS: mean= {AmQRS}, average amplitude of fQRS: mean= {AfQRS}')
print(f'ratio maternal/fetal : {AmQRS/AfQRS}, ratio fetal/maternal: {AfQRS/AmQRS} ')

average amplitude of mQRS: mean= 103.25267355520559, average amplitude of fQRS: mean= 20.49577544888795
ratio maternal/fetal : 5.03775394166937, ratio fetal/maternal: 0.19850115975863403

Question 4¶

How many of the principal components of the first PCA clearly show a maternal ECG signal? Which ones?

In [15]:
print(pca1.explained_variance_ratio_)
fig = go.Figure()

fig.add_trace(go.Scatter(x=[1,2,3,4], y=pca1.explained_variance_ratio_, mode='lines'))

fig.update_xaxes(title='components')
fig.update_yaxes(title='Variance ratio')
fig.update_layout(title='Scree plot of the maternal ECG signal')

fig.show()

[0.65448309 0.32978623 0.01279188 0.0029388 ]

We can see on the plot above that we are seeing the first two components of the first PCA.

first PCA first PCA

Question 5¶

How many of the principal components of the second PCA clearly show a fetal ECG signal? Which ones?

In [23]:
print(pca2.explained_variance_ratio_)
fig = go.Figure()

fig.add_trace(go.Scatter(x=[1,2,3,4], y=pca2.explained_variance_ratio_, mode='lines'))

fig.update_xaxes(title='components')
fig.update_yaxes(title='Variance ratio')
fig.update_layout(title='Scree plot of the fetal ECG signal')

fig.show()

[0.77164836 0.12371719 0.0672538  0.03738065]

We can see on the plot above that we are seeing the first component of the second PCA of the fetal ECG signal.

second PCA second PCA

Question 6¶

Not all of the fetal QRS peaks seem to be detected properly. Do you have an explanation why this happens and under which circumstances? Is it a problem of the fQRS detector, the mQRS cancellation or of another block of the algorithm?

One of the possible reasons that some peaks of the fetal QRS are not detected can be that they have a too closed proximity with the peaks of the mother QRS. Therefore, it could be that the cancellation is not perfect and we still have some residues of mQRS in the fQRS signal.

comp